We consider the Hamiltonian
$$ H(r, \theta, p_r, p_\theta) = v p_\theta + \Vert p \Vert_{g} $$
where $v$ is a constant, $p = (p_r, p_\theta)$, and $\Vert \cdot \Vert_{g}$ is the norm induced by the metric
$$ g = \mathrm{d}r^2 + m_\lambda^2(r)\, \mathrm{d}\theta^2, \quad m_\lambda^2(r) = \frac{\sin^2 r}{1 - \lambda \sin^2 r} $$
with $\lambda = 4/5$.
Along the geodesics, we have $H+p^0 = 0$. The parameter $p^0$ is constant equal to $-1$ (hyperbolic), $0$ (abnormal) or $1$ (elliptic).
Remark. We can parameterize the geodesics by the norm of the initial convector, setting $\Vert{p_0}\Vert_g = 1$. This amounts to parameterize by the initial angle $\alpha_0$: $$ p_r = \sin \alpha_0, \quad p_\theta = m_\lambda(r) \cos \alpha_0. $$ In that case, the hyperbolic geodeics are given by $$ p_\theta\, v + 1 = v\, m_\lambda(r) \cos \alpha_0 + 1 > 0. $$
import sys
sys.path.insert(1, '../')
#
import math
import numpy as np # scientific computing tools
import wrappers # for compilation of Fortran Hamiltonian codes
import geometry2d.database
import geometry2d.conjugate
import geometry2d.errors
import geometry2d.geodesic
import geometry2d.plottings
import geometry2d.problem
import geometry2d.splitting
import geometry2d.utils
import geometry2d.wavefront
import geometry2d.sphere
import nutopy as nt
import scipy
import matplotlib.pyplot as plt # for plots
%matplotlib inline
#
color_hyperbolic_2D = 'red'
color_hyperbolic_3D = (1.0, 0.2, 0.2)
color_elliptic_2D = 'blue'
color_elliptic_3D = (0.2, 0.2, 1.0)
color_abnormal_2D = (0.2, 0.8, 0.2) #'green'
color_abnormal_3D = (0.2, 1.0, 0.2)
color_cut_locus = 'black'
color_conjugate = 'magenta'
color_wavefront = 'orange'
color_sphere = 'orange'
#
color_strong_domain_2d = 'lightblue'
color_strong_domain_3d = 'black'
#
PLANE = geometry2d.plottings.Coords.PLANE
SPHERE = geometry2d.plottings.Coords.SPHERE
# print Fortran code of the Hamiltonian: print hfun.f90 file
with open('hfun.f90', 'r') as f:
print(f.read())
! Kepler
subroutine hfun(x, p, v, h)
double precision, intent(in) :: x(2), p(2), v
double precision, intent(out) :: h
! local variables
double precision :: r, th, pr, pth
double precision :: l, m
r = x(1)
th = x(2)
pr = p(1)
pth = p(2)
l = 4d0/5d0
m = sqrt(sin(r)**2/(1d0-l*sin(r)**2))
h = v*pth + sqrt(pr**2 + (pth/m)**2)
end subroutine hfun
def m2(r, λ):
return np.sin(r)**2/(1.0-λ*np.sin(r)**2)
λ_Kepler = 4.0/5.0
def g2_Kepler(r):
return m2(r, λ_Kepler)
def initialize(v, r0):
# Parameters
θ0 = 0.0 # initial latitude
t0 = 0.0 # initial time
name = 'kepler_strong' # name of the problem
# Initialize data
data_file = 'data_kepler_strong.json'
restart = False # restart or not the computations
data = geometry2d.database.Data({'name': name,
't0': t0,
'r0': r0,
'θ0': θ0,
'v': v}, data_file, restart)
# Initial point
q0 = np.array([r0, θ0])
# Hamiltonian and derivatives up to order 3
H = wrappers.hamiltonian(v, compile=False, display=False)
# The Riemannian metric associated to the Zermelo problem
# g = g1 dr^2 + g2 dθ^2
def g(q):
#λ = 4.0/5.0
r = q[0]
g1 = 1.0
g2 = g2_Kepler(r)
return g1, g2
# problem
prob = geometry2d.problem.GeometryProblem2D(name, H, g, t0, q0, data, steps_for_geodesics=200)
return prob
# we parameterize ||p(0)|| = 1 by the angle α0
# the norm being given by the metric g: p(0) = (sin(α0)*g1, cos(α0)*g2) = (sin(α0), cos(α0)*g2), since g1 = 1
# get the values of α0 for which h = ptheta v + 1 > 0
#
# first get the values of α for which h = ptheta v + 1 = 0
# since ptheta = m(r) cos(α0) we wil need to check if v m(r) > 1
# m(r) = sqrt(g2_Kepler(r))
def K_fun(r, v):
return v*np.sqrt(g2_Kepler(r))
def pθ_fun(r, α):
return np.cos(α)*np.sqrt(g2_Kepler(r))
def h_fun(r, α, v):
return pθ_fun(r, α)*v + 1.0
def α0_fun(r, v):
return np.arccos(-1.0/K_fun(r, v))
# get the values of α0 for which h > 0
def α_hyperbolique_1(r, v):
k = K_fun(r, v)
if abs(k) > 1.0:
α = α0_fun(r, v) # between 0 and pi
if α < np.pi/2.0:
return α
else:
return 2.0*np.pi - α
else:
return 0.0
def α_hyperbolique_2(r, v):
k = K_fun(r, v)
if abs(k) > 1.0:
α = α0_fun(r, v) # between 0 and pi
if α < np.pi/2.0:
return 2.0*np.pi - α
else:
return α + 2.0*np.pi
else:
return 2.0*np.pi
def α_hyperbolique_span(r, v, *, N=20):
α1 = α_hyperbolique_1(r, v)
print('α1 = ', α1)
α2 = α_hyperbolique_2(r, v)
print('α2 = ', α2)
return np.linspace(α1, α2, N)
# get the values of α0 for which h < 0
def α_elliptique_1(r, v):
k = K_fun(r, v)
if abs(k) > 1.0:
α = α0_fun(r, v) # between 0 and pi
if α < np.pi/2.0:
return 2.0*np.pi - α
else:
return α
else:
# raise an error: there is no elliptic domain
raise geometry2d.errors.ArgumentValueError('There is no elliptic domain for this value of v')
def α_elliptique_2(r, v):
k = K_fun(r, v)
if abs(k) > 1.0:
α = α0_fun(r, v) # between 0 and pi
if α < np.pi/2.0:
return α + 2.0*np.pi
else:
return 2.0*np.pi - α
else:
# raise an error: there is no elliptic domain
raise geometry2d.errors.ArgumentValueError('There is no elliptic domain for this value of v')
def α_elliptique_span(r, v, *, N=20):
α1 = α_elliptique_1(r, v)
print('α1 = ', α1)
α2 = α_elliptique_2(r, v)
print('α2 = ', α2)
return np.linspace(α1, α2, N)
def α_abnormals(r, v):
k = K_fun(r, v)
if abs(k) >= 1.0:
α = α0_fun(r, v)
if α < np.pi/2.0:
return np.array([α, -α])
else:
return np.array([α, 2.0*np.pi-α])
else:
# raise an error: there is no abnormals geodesics
raise geometry2d.errors.ArgumentValueError('There is no abnormals geodesics for this value of v')
def domain_strong_current(v):
r1 = math.asin(np.sqrt(1.0/(λ_Kepler+v**2)))
r2 = np.pi - r1
return r1, r2
def update_2d_plot(fig):
ax = fig.axes[0]
# set xlims
ax.set_xlim(-np.pi/3, 2*np.pi + np.pi/3)
# add vertical lines
ax.axvline(2*np.pi, color='k', linewidth=0.5, linestyle="dashed", zorder=geometry2d.plottings.z_order_axes)
ax.axvline(1*np.pi, color='k', linewidth=0.5, linestyle="dashed", zorder=geometry2d.plottings.z_order_axes)
# set ylabel to "r"
ax.set_ylabel(r'$r$', rotation=0, labelpad=10)
# set yticks: 0 => \pi/2, pi/2 => \pi, -pi/2 => 0
ax.set_yticks([-np.pi/2, 0.0, np.pi/2])
ax.set_yticklabels([r'$0$', r'$\frac{\pi}{2}$', r'$\pi$'])
# set xticks: 0, pi and 2pi
ax.set_xticks([0.0, np.pi, 2*np.pi])
ax.set_xticklabels([r'$0$', r'$\pi$', r'$2\pi$'])
return fig
v = 0.8
r0 = np.pi/2.0
prob = initialize(v, r0)
geodesic = geometry2d.geodesic.Geodesic(prob)
r1_strong, r2_strong = domain_strong_current(v) # Domain of strong current
v_Kepler = v
αspan_hyperbolique = α_hyperbolique_span(r0, v, N=21)
tf = geodesic.return_to_equator
#
fig2d = geodesic.plot(alphas=αspan_hyperbolique, tf=tf, length=1.0, view=PLANE, figsize=(5,5), color=color_hyperbolic_2D)
fig2d = update_2d_plot(fig2d)
fig2d = geometry2d.utils.plot_2d_domain(fig2d, r1_strong, r2_strong, color_strong_domain_2d)
# cameras is a list of dictionaries containing azimuth and elevetion angles
cameras=[{'azimuth': 140, 'elevation': -10}]
figs3d = []
for cam in cameras:
azimuth = cam['azimuth']
elevation = cam['elevation']
fig3d = geodesic.plot(alphas=αspan_hyperbolique, tf=tf, length=1.0, view=SPHERE,
azimuth=azimuth, elevation=elevation, figsize=(3,3), color=color_hyperbolic_3D)
fig3d = geometry2d.utils.plot_3d_domain(fig3d, r1_strong, r2_strong, color_strong_domain_3d, azimuth=azimuth, elevation=elevation)
figs3d.append(fig3d)
α1 = 4.119189204235061 α2 = 8.447181410124111
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
LW_GEO = 0.8
α1a, α2a = α_abnormals(r0, v)
print('α1a = ', α1a)
print('α2a = ', α2a)
α1a = 2.1639961029445254 α2a = 4.119189204235061
tf = geodesic.return_to_equator
#
fig2d = geodesic.plot(alphas=[α1a, α2a], tf=tf, linewidth=LW_GEO, view=PLANE, figsize=(5,5), color=color_abnormal_2D)
fig2d = update_2d_plot(fig2d)
fig2d = geometry2d.utils.plot_2d_domain(fig2d, r1_strong, r2_strong, color_strong_domain_2d)
# cameras is a list of dictionaries containing azimuth and elevetion angles
cameras=[{'azimuth': 140, 'elevation': -10}]
figs3d = []
for cam in cameras:
azimuth = cam['azimuth']
elevation = cam['elevation']
fig3d = geodesic.plot(alphas=[α1a, α2a], tf=tf, linewidth=LW_GEO, view=SPHERE,
azimuth=azimuth, elevation=elevation, figsize=(3,3), color=color_abnormal_3D)
fig3d = geometry2d.utils.plot_3d_domain(fig3d, r1_strong, r2_strong, color_strong_domain_3d, azimuth=azimuth, elevation=elevation)
figs3d.append(fig3d)
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
# the separatrices between hyperbolic with loops and hyperbolic without loops is given by pθ = 0
# without loops:
αs_no_loops = np.linspace(0.0, np.pi/2, 5)
αs_no_loops = αs_no_loops[[1,3]]
# with loops
αe1 = α_elliptique_1(r0, v)
αs_loops = np.linspace(np.pi/2, αe1, 7)
αs_loops = αs_loops[[1,4]]
# concatenate
αs = np.concatenate((αs_no_loops, αs_loops))
#
tf = geodesic.return_to_equator
#
fig2d = geodesic.plot(alphas=αs, tf=tf, linewidth=LW_GEO, view=PLANE, color=color_hyperbolic_2D, figure=fig2d)
#
for fig3d in figs3d:
geodesic.plot(alphas=αs, tf=tf, linewidth=LW_GEO, view=SPHERE, color=color_hyperbolic_3D, figure=fig3d)
fig2d
figs3d[0]
# get the values of α0 for which the hyperbolic geodesics return to the initial point
def myfun_return_to_initial(y):
t, α0 = y
q = prob.geodesic(t, α0)
return q[0] - q[-1]
t_loop_guess = 3.0
# α0_loop_guess = middle of pi/2 and abnormal given by α1a
α0_loop_guess = 1.7 #(np.pi/2 + α1a)/2.0
y_guess = np.array([t_loop_guess, α0_loop_guess])
sol = nt.nle.solve(myfun_return_to_initial, y_guess)
t_loop_return_to_initial = sol.x[0]
α_loop_return_to_initial_1 = sol.x[1]
α_loop_return_to_initial_2 = 2*np.pi - α_loop_return_to_initial_1
Calls |f(x)| |x|
1 5.651442610204609e-02 3.448187929913334e+00
2 3.055349141119294e-05 3.483630790700333e+00
3 3.106629222756182e-07 3.483618968834663e+00
4 7.495891880478653e-11 3.483618817869363e+00
5 1.404519055154983e-14 3.483618817899530e+00
Results of the nle solver method:
xsol = [3.03970007 1.70171195]
f(xsol) = [-6.88338275e-15 -1.22428109e-14]
nfev = 5
njev = 1
status = 1
success = True
Successfully completed: relative error between two consecutive iterates is at most TolX.
#
αe1 = α_elliptique_1(r0, v)
αs = np.linspace(αe1, np.pi, 4)
αs = αs[1:-1]
#
tf = geodesic.return_to_equator
#
fig2d = geodesic.plot(alphas=αs, tf=tf, linewidth=LW_GEO, view=PLANE, color=color_elliptic_2D, figure=fig2d)
#
for fig3d in figs3d:
geodesic.plot(alphas=αs, tf=tf, linewidth=LW_GEO, view=SPHERE, color=color_elliptic_3D, figure=fig3d)
fig2d.savefig('kepler_geodesic_flow_2D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
fig2d
figs3d[0].savefig('kepler_geodesic_flow_3D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
figs3d[0]
#
gap = 1e-3
# hyperbolic whitout loops
αh_no_loops_1 = -np.pi/2 + gap
αh_no_loops_2 = np.pi/2 - gap
# hyperbolic with loops
αh_loops_1 = np.pi/2 + gap
αh_loops_2 = α_elliptique_1(r0, v) - gap
αh_loops_3 = α_elliptique_2(r0, v) + gap
αh_loops_4 = 3.0*np.pi/2 - gap
# elliptic
αe_1 = α_elliptique_1(r0, v) + gap
αe_2 = α_elliptique_2(r0, v) - gap
conjugate = geometry2d.conjugate.Conjugate(prob)
conj_h_no_loops = conjugate.compute([αh_no_loops_1, αh_no_loops_2], label='hyperbolic_without_loops')
conj_h_loop_north = conjugate.compute([αh_loops_1, αh_loops_2], label='hyperbolic_with_loops_north')
conj_h_loop_south = conjugate.compute([αh_loops_3, αh_loops_4], label='hyperbolic_with_loops_south')
conj_elliptic = conjugate.compute([αe_1, αe_2], label='elliptic')
LW_CONJ = 1.5
# PLANE
fig2d = conjugate.plot(['hyperbolic_without_loops', 'hyperbolic_with_loops_north', 'hyperbolic_with_loops_south'],
view=PLANE, color=color_conjugate, figure=fig2d, linestyle='solid', linewidth=LW_CONJ)
fig2d = conjugate.plot(['elliptic'], view=PLANE, color=color_conjugate, figure=fig2d, linestyle='dashed', linewidth=LW_CONJ)
# SPHERE
for fig3d in figs3d:
conjugate.plot(['hyperbolic_without_loops', 'hyperbolic_with_loops_north', 'hyperbolic_with_loops_south'],
view=SPHERE, color=color_conjugate, figure=fig3d, linestyle='solid', linewidth=LW_CONJ)
conjugate.plot(['elliptic'], view=SPHERE, color=color_conjugate, figure=fig3d, linestyle='dashed', linewidth=LW_CONJ)
fig2d.savefig('kepler_geodesic_flow_2D_with_conjugate_locus.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
fig2d
figs3d[0].savefig('kepler_geodesic_flow_3D_with_conjugate_locus.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
figs3d[0]
wavefront = geometry2d.wavefront.WaveFront(prob)
tfs_wavefronts = np.arange(0.1, 3.5, 0.1)
if False:
for tf in tfs_wavefronts:
wavefront.compute(tf)
tfs = [2.5, 3.1, 3.4]
for i in range(len(tfs)):
tf = tfs[i]
if i == 0:
fig2d_wavefront = wavefront.plot(tf, view=PLANE, color=color_wavefront, figsize=(5,5))
fig2d_wavefront = update_2d_plot(fig2d_wavefront)
fig2d_wavefront = geometry2d.utils.plot_2d_domain(fig2d_wavefront, r1_strong, r2_strong, color_strong_domain_2d)
figs3d_wavefront = []
for camera in cameras:
azimuth = camera['azimuth']
elevation = camera['elevation']
fig3d_wavefront = wavefront.plot(tf, view=SPHERE, color=color_wavefront,
azimuth=azimuth, elevation=elevation, figsize=(3,3))
fig3d_wavefront = geometry2d.utils.plot_3d_domain(fig3d_wavefront, r1_strong, r2_strong, color_strong_domain_3d, azimuth=azimuth, elevation=elevation)
figs3d_wavefront.append(fig3d_wavefront)
else:
wavefront.plot(tf, view=PLANE, color=color_wavefront, figure=fig2d_wavefront)
for fig3d_wavefront in figs3d_wavefront:
wavefront.plot(tf, view=SPHERE, color=color_wavefront, figure=fig3d_wavefront)
# add conjugate locus
LW_CONJ = 1.5
# PLANE
fig2d_wavefront = conjugate.plot(['hyperbolic_without_loops', 'hyperbolic_with_loops_north', 'hyperbolic_with_loops_south'],
view=PLANE, color=color_conjugate, figure=fig2d_wavefront, linestyle='solid', linewidth=LW_CONJ)
fig2d_wavefront = conjugate.plot(['elliptic'], view=PLANE, color=color_conjugate, figure=fig2d_wavefront, linestyle='dashed', linewidth=LW_CONJ)
# SPHERE
for fig3d in figs3d_wavefront:
conjugate.plot(['hyperbolic_without_loops', 'hyperbolic_with_loops_north', 'hyperbolic_with_loops_south'],
view=SPHERE, color=color_conjugate, figure=fig3d, linestyle='solid', linewidth=LW_CONJ)
conjugate.plot(['elliptic'], view=SPHERE, color=color_conjugate, figure=fig3d, linestyle='dashed', linewidth=LW_CONJ)
<Figure size 640x480 with 0 Axes>
fig2d_wavefront.savefig('kepler_wavefronts_2D_with_conjugate_locus.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
fig2d_wavefront
figs3d_wavefront[0].savefig('kepler_wavefronts_2D_with_conjugate_locus.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
figs3d_wavefront[0]
splitting = geometry2d.splitting.Splitting(prob)
# no loops
intersections = wavefront.self_intersections(2.0)
print(intersections['left'])
print(intersections['right'])
intersections_no_loops = intersections['right']
αspan = [0+gap, αh_no_loops_2]
print(αspan)
cut_no_loops = splitting.compute(intersections_no_loops, αspan, label='no_loops')
t_cut_no_loops = splitting.splitting_time(cut_no_loops)
(2.0, 2.2214256373524743, 4.0617516058063154, array([1.57079826, 0.62535495])) (2.0, -0.9201320258729991, 0.9201777291167127, array([1.57080738, 2.57464891])) [0.001, 1.5697963267948967]
# loops: α in [α_loop_return_to_initial_1, pi/2]
intersections = wavefront.self_intersections(3.1)
print(intersections['left'])
print(intersections['right'])
intersections_loops = intersections['left']
t = intersections_loops[0]
α1 = intersections_loops[1]
α2 = intersections_loops[2]
q = intersections_loops[3]
intersections_loops = (t, α2, α1, q)
αspan = [α_loop_return_to_initial_1, np.pi/2+gap]
print(αspan)
cut_loops = splitting.compute(intersections_loops, αspan, label='loops')
t_cut_loops = splitting.splitting_time(cut_loops)
(3.1, 1.6530290659031264, 4.63015253889408, array([ 1.57081293, -0.20609652])) (3.1, -1.4886017740158815, 1.4885861368789068, array([1.5707231 , 5.16628826])) [1.7017119516020776, 1.5717963267948964]
print(α_loop_return_to_initial_2)
cut_loops.alphas
4.581473355577509
array([[4.58147336, 1.70171195],
[4.58240506, 1.70078024],
[4.58962422, 1.69356109],
[4.60121676, 1.68196855],
[4.61711397, 1.66607133],
[4.63632337, 1.64686193],
[4.6582413 , 1.62494401],
[4.68249347, 1.60069184],
[4.70907165, 1.57411366],
[4.71138898, 1.57179633]])
LW_CUT = 1.8
# plot cut locus
fig2d_synthesis = splitting.plot(['no_loops', 'loops'], linewidth=LW_CUT, view=PLANE, color=color_cut_locus,
figsize=(5,5))
fig2d_synthesis = update_2d_plot(fig2d_synthesis)
fig2d_synthesis = geometry2d.utils.plot_2d_domain(fig2d_synthesis, r1_strong, r2_strong, color_strong_domain_2d)
cameras=[{'azimuth': 140, 'elevation': -10}]
figs3d_synthesis = []
for cam in cameras:
azimuth = cam['azimuth']
elevation = cam['elevation']
fig3d_synthesis = splitting.plot(['no_loops', 'loops'], linewidth=LW_CUT, view=SPHERE, color=color_cut_locus,
azimuth=azimuth, elevation=elevation, figsize=(3,3))
fig3d_synthesis = geometry2d.utils.plot_3d_domain(fig3d_synthesis, r1_strong, r2_strong, color_strong_domain_3d, azimuth=azimuth, elevation=elevation)
figs3d_synthesis.append(fig3d_synthesis)
<Figure size 640x480 with 0 Axes>
LW_CONJ = 1.5
# add hyperbolic conjugate points
fig2d_synthesis = conjugate.plot(['hyperbolic_with_loops_north', 'hyperbolic_with_loops_south',
'hyperbolic_without_loops'], linewidth=LW_CONJ, view=PLANE,
color=color_conjugate, figure=fig2d_synthesis, linestyle='solid')
for fig3d_synthesis in figs3d_synthesis:
conjugate.plot(['hyperbolic_with_loops_north', 'hyperbolic_with_loops_south',
'hyperbolic_without_loops'], linewidth=LW_CONJ, view=SPHERE,
color=color_conjugate, figure=fig3d_synthesis, linestyle='solid')
# add abnormal geodesics until the cusp point defined by the boundary of the domain of strong current
α1a, α2a = α_abnormals(r0, v)
#
print('α1a = ', α1a)
print('α2a = ', α2a)
print('r1_strong = ', r1_strong)
print('r2_strong = ', r2_strong)
# the abnormal associated to α1a meets the boundary r = r2_strong
def myfun(t):
return prob.geodesic(t, α1a)[-1][0] - (r2_strong - 1e-5)
opt = nt.nle.Options(Display='off')
t_cusp_α1a = nt.nle.solve(myfun, 1.0, options=opt).x
# the abnormal associated to α2a meets the boundary r = r1_strong
def myfun(t):
return prob.geodesic(t, α2a)[-1][0] - (r1_strong + 1e-5)
t_cusp_α2a = nt.nle.solve(myfun, 1.0, options=opt).x
# t_cusp_abnormal
def t_cusp_abnormal(α):
if (α % 2*np.pi) == (α1a % 2*np.pi):
return t_cusp_α1a
elif (α % 2*np.pi) == (α2a % 2*np.pi):
return t_cusp_α2a
else:
raise geometry2d.errors.ArgumentValueError('α must be equal to α1a or α2a')
# plot
fig2d_synthesis = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=LW_CUT,
color=color_abnormal_2D, figure=fig2d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
fig2d_synthesis = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=LW_CUT,
color=color_cut_locus, linestyle='dotted', figure=fig2d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
for fig3d_synthesis in figs3d_synthesis:
geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=LW_CUT,
color=color_abnormal_3D, figure=fig3d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=LW_CUT,
color=color_cut_locus, linestyle='dotted', figure=fig3d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
α1a = 2.1639961029445254 α2a = 4.119189204235061 r1_strong = 0.9851107833377455 r2_strong = 2.1564818702520476
fig2d_synthesis.savefig('kepler_conj_cut_2D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
fig2d_synthesis
figs3d_synthesis[0].savefig('kepler_conj_cut_3D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
figs3d_synthesis[0]
# plot cut locus
fig2d_synthesis = splitting.plot(['no_loops', 'loops'], linewidth=LW_CUT, view=PLANE, color=color_cut_locus,
figsize=(5,5))
fig2d_synthesis = update_2d_plot(fig2d_synthesis)
fig2d_synthesis = geometry2d.utils.plot_2d_domain(fig2d_synthesis, r1_strong, r2_strong, color_strong_domain_2d)
cameras=[{'azimuth': 140, 'elevation': -10}]
figs3d_synthesis = []
for cam in cameras:
azimuth = cam['azimuth']
elevation = cam['elevation']
fig3d_synthesis = splitting.plot(['no_loops', 'loops'], linewidth=LW_CUT, view=SPHERE, color=color_cut_locus,
azimuth=azimuth, elevation=elevation, figsize=(3,3))
fig3d_synthesis = geometry2d.utils.plot_3d_domain(fig3d_synthesis, r1_strong, r2_strong, color_strong_domain_3d, azimuth=azimuth, elevation=elevation)
figs3d_synthesis.append(fig3d_synthesis)
# add abnormal geodesics until the cusp point defined by the boundary of the domain of strong current
α1a, α2a = α_abnormals(r0, v)
#
print('α1a = ', α1a)
print('α2a = ', α2a)
print('r1_strong = ', r1_strong)
print('r2_strong = ', r2_strong)
# the abnormal associated to α1a meets the boundary r = r2_strong
def myfun(t):
return prob.geodesic(t, α1a)[-1][0] - (r2_strong - 1e-5)
opt = nt.nle.Options(Display='off')
t_cusp_α1a = nt.nle.solve(myfun, 1.0, options=opt).x
# the abnormal associated to α2a meets the boundary r = r1_strong
def myfun(t):
return prob.geodesic(t, α2a)[-1][0] - (r1_strong + 1e-5)
t_cusp_α2a = nt.nle.solve(myfun, 1.0, options=opt).x
# t_cusp_abnormal
def t_cusp_abnormal(α):
if (α % 2*np.pi) == (α1a % 2*np.pi):
return t_cusp_α1a
elif (α % 2*np.pi) == (α2a % 2*np.pi):
return t_cusp_α2a
else:
raise geometry2d.errors.ArgumentValueError('α must be equal to α1a or α2a')
# plot
fig2d_synthesis = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=LW_CUT,
color=color_abnormal_2D, figure=fig2d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
fig2d_synthesis = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=LW_CUT,
color=color_cut_locus, linestyle='dotted', figure=fig2d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
for fig3d_synthesis in figs3d_synthesis:
geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=LW_CUT,
color=color_abnormal_3D, figure=fig3d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=LW_CUT,
color=color_cut_locus, linestyle='dotted', figure=fig3d_synthesis,
zorder=geometry2d.plottings.z_order_splitting)
α1a = 2.1639961029445254 α2a = 4.119189204235061 r1_strong = 0.9851107833377455 r2_strong = 2.1564818702520476
<Figure size 640x480 with 0 Axes>
restart_path_loop_1 = False
if prob.data.contains('path_loop_1') and not restart_path_loop_1:
path_loop_1 = prob.data.get('path_loop_1')
yout_loop_1 = np.array(path_loop_1['yout'])
αout_loop_1 = np.array(path_loop_1['αout'])
print('path_loop_1 loaded from data')
else:
# add hyperbolic geodesics with loops until the cut point defined by the intersection with the abnormals
# for α in [pi/2, αh_loops_1] the intersection is with the abnormal α1a
def myfun_itersection_1(y, α):
th, ta, r, θ = y
q = np.array([r, θ])
eq = np.zeros(4)
qf = prob.geodesic(th, α[0])[-1]
eq[0:2] = q - qf
eq[2:4] = q - prob.geodesic(ta, α1a)[-1]
return eq
# we start with α close to α1a since we know a good approximation of the intersection
α0_ = α1a - 1e-3
th = t_cusp_α1a
ta = t_cusp_α1a
q = prob.geodesic(th, α0_)[-1]
r = q[0]
θ = q[1]
y_guess = np.array([th, ta, r, θ])
sol = nt.nle.solve(lambda y : myfun_itersection_1(y, [α0_]), y_guess)
y0 = sol.x
print('th = ', y0[0])
print('ta = ', y0[1])
print('t_cusp_α1a = ', t_cusp_α1a)
sol = nt.path.solve(myfun_itersection_1, y0, [α0_], [α_loop_return_to_initial_1])
yout_loop_1, αout_loop_1 = sol.xout, sol.parsout
path_loop_1 = {'yout': yout_loop_1.tolist(), 'αout': αout_loop_1.tolist()}
prob.data.update({'path_loop_1': path_loop_1})
print('path_loop_1 saved in data')
path_loop_1 loaded from data
restart_path_loop_2 = False
if prob.data.contains('path_loop_2') and not restart_path_loop_2:
path_loop_2 = prob.data.get('path_loop_2')
yout_loop_2 = np.array(path_loop_2['yout'])
αout_loop_2 = np.array(path_loop_2['αout'])
print('path_loop_2 loaded from data')
else:
# add hyperbolic geodesics with loops until the cut point defined by the intersection with the abnormals
# for α in [pi/2, αh_loops_1] the intersection is with the abnormal α1a
def myfun_itersection_2(y, α):
th, ta, r, θ = y
q = np.array([r, θ])
eq = np.zeros(4)
qf = prob.geodesic(th, α[0])[-1]
eq[0:2] = q - qf
eq[2:4] = q - prob.geodesic(ta, α2a)[-1]
return eq
# we start with α close to α1a since we know a good approximation of the intersection
α0_ = α2a + 1e-3
th = t_cusp_α2a
ta = t_cusp_α2a
q = prob.geodesic(th, α0_)[-1]
r = q[0]
θ = q[1]
y_guess = np.array([th, ta, r, θ])
sol = nt.nle.solve(lambda y : myfun_itersection_2(y, [α0_]), y_guess)
y0 = sol.x
print('th = ', y0[0])
print('ta = ', y0[1])
print('t_cusp_α2a = ', t_cusp_α2a)
sol = nt.path.solve(myfun_itersection_2, y0, [α0_], [α_loop_return_to_initial_2])
yout_loop_2, αout_loop_2 = sol.xout, sol.parsout
path_loop_2 = {'yout': yout_loop_2.tolist(), 'αout': αout_loop_2.tolist()}
prob.data.update({'path_loop_2': path_loop_2})
print('path_loop_2 saved in data')
path_loop_2 loaded from data
LW_GEO = 0.8
# add hyperbolic geodesics with loops until the cut point defined by the intersection with the abnormals
def t_meet_abnormal(α):
α = α % (2.0*np.pi)
if α <= (α1a % (2.0*np.pi)) and α >= (α_loop_return_to_initial_1 % (2.0*np.pi)):
times = yout_loop_1[:, 0]
alphas = αout_loop_1[:, 0] % (2*np.pi)
elif α >= (α2a % (2.0*np.pi)) and α <= (α_loop_return_to_initial_2 % (2.0*np.pi)):
times = yout_loop_2[:, 0]
alphas = αout_loop_2[:, 0] % (2*np.pi)
else:
raise geometry2d.errors.ArgumentValueError('α must be in [α_loop_return_to_initial_1, α1a] U [α2a, α_loop_return_to_initial_2]')
t_meet = scipy.interpolate.interp1d(alphas, times, kind='cubic', fill_value="extrapolate")(α)
t_equator = geodesic.first_return_to_equator(α)
return min(t_meet, t_equator)
# plot
αspan1 = np.linspace(α_loop_return_to_initial_1+1e-3, α1a-1e-3, 5)
αspan1 = αspan1[0:-1]
αspan2 = np.linspace(α_loop_return_to_initial_2-1e-3, α2a+1e-3, 5)
αspan2 = αspan2[0:-1]
αspan = np.concatenate((αspan1, αspan2))
fig2d_synthesis = geodesic.plot(alphas=αspan, tf=t_meet_abnormal, view=PLANE, linewidth=LW_GEO,
color=color_hyperbolic_2D, figure=fig2d_synthesis)
for fig3d_synthesis in figs3d_synthesis:
geodesic.plot(alphas=αspan, tf=t_meet_abnormal, view=SPHERE, linewidth=LW_GEO,
color=color_hyperbolic_3D, figure=fig3d_synthesis)
fig2d_synthesis
# add hyperbolic geodesics without loops until the cut locus
αspan = np.linspace(αh_no_loops_1, αh_no_loops_2, 11)
αspan = np.concatenate(( αspan, [(αspan[-1]+αspan[-2])/2.0], [(αspan[0]+αspan[1])/2.0] ))
fig2d_synthesis = geodesic.plot(alphas=αspan, tf=t_cut_no_loops, view=PLANE, linewidth=LW_GEO,
color=color_hyperbolic_2D, figure=fig2d_synthesis)
for fig3d_synthesis in figs3d_synthesis:
geodesic.plot(alphas=αspan, tf=t_cut_no_loops, view=SPHERE, linewidth=LW_GEO,
color=color_hyperbolic_3D, figure=fig3d_synthesis)
fig2d_synthesis
# add hyperbolic geodesics with loops but not intersection with abnormals until the cut locus
αspan1 = np.linspace(np.pi/2+1e-3, α_loop_return_to_initial_1, 3)
αspan2 = np.linspace(α_loop_return_to_initial_2+1e-3, 3.0*np.pi/2-1e-3, 3)
αspan = np.concatenate(( αspan1, αspan2 ))
fig2d_synthesis = geodesic.plot(alphas=αspan, tf=t_cut_loops, view=PLANE, linewidth=LW_GEO,
color=color_hyperbolic_2D, figure=fig2d_synthesis)
for fig3d_synthesis in figs3d_synthesis:
geodesic.plot(alphas=αspan, tf=t_cut_loops, view=SPHERE, linewidth=LW_GEO,
color=color_hyperbolic_3D, figure=fig3d_synthesis)
fig2d_synthesis.savefig('kepler_synthesis_2D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
fig2d_synthesis
figs3d_synthesis[0].savefig('kepler_synthesis_3D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
figs3d_synthesis[0]

def t_cut_final(α):
α = α % (2.0*np.pi)
α1a, α2a = α_abnormals(r0, v)
α1a = α1a % (2.0*np.pi)
α2a = α2a % (2.0*np.pi)
α1_loop = np.pi/2
α2_loop = 3.0*np.pi/2
α1_loop_return = α_loop_return_to_initial_1 % (2.0*np.pi)
α2_loop_return = α_loop_return_to_initial_2 % (2.0*np.pi)
if ( 0.0 <= α <= α1_loop ) or ( α2_loop <= α <= 2.0*np.pi ):
return t_cut_no_loops(α)
elif ( α1_loop <= α <= α1_loop_return ) or ( α2_loop_return <= α <= α2_loop ): # meet abnormal
return t_cut_loops(α)
elif ( α1_loop_return <= α < α1a ) or ( α2a < α <= α2_loop_return ): # loops
return t_meet_abnormal(α)
elif α1a == α or α2a == α:
return t_cusp_abnormal(α)
else:
return 0.0
sphere = geometry2d.sphere.Sphere(wavefront, t_cut_final)
if False:
for tf in tfs_wavefronts:
print(tf)
sphere.compute(tf)
LW_SPHERE_2D = 0.8
LW_SPHERE_3D = 1.2
tfs = tfs_wavefronts[0:-1:2]
for i in range(len(tfs)):
tf = tfs[i]
if i == 0:
fig2d_sphere = sphere.plot(tf, view=PLANE, color=color_sphere, figsize=(5,5), linewidth=LW_SPHERE_2D)
fig2d_sphere = update_2d_plot(fig2d_sphere)
fig2d_sphere = geometry2d.utils.plot_2d_domain(fig2d_sphere, r1_strong, r2_strong, color_strong_domain_2d)
figs3d_sphere = []
for camera in cameras:
azimuth = camera['azimuth']
elevation = camera['elevation']
fig3d_sphere = sphere.plot(tf, view=SPHERE, color=color_sphere, linewidth=LW_SPHERE_3D,
azimuth=azimuth, elevation=elevation, figsize=(3,3))
fig3d_sphere = geometry2d.utils.plot_3d_domain(fig3d_sphere, r1_strong, r2_strong, color_strong_domain_3d, azimuth=azimuth, elevation=elevation)
figs3d_sphere.append(fig3d_sphere)
else:
sphere.plot(tf, view=PLANE, color=color_sphere, figure=fig2d_sphere, linewidth=LW_SPHERE_2D)
for fig3d_sphere in figs3d_sphere:
sphere.plot(tf, view=SPHERE, color=color_sphere, figure=fig3d_sphere, linewidth=LW_SPHERE_3D)
<Figure size 640x480 with 0 Axes>
LW_CUT = 1.8
# plot cut locus
fig2d_sphere = splitting.plot(['no_loops', 'loops'], linewidth=LW_CUT, view=PLANE, figure=fig2d_sphere,
color=color_cut_locus)
for fig3d_sphere in figs3d_sphere:
fig3d_sphere = splitting.plot(['no_loops', 'loops'], linewidth=LW_CUT, view=SPHERE, figure=fig3d_sphere,
color=color_cut_locus)
# # add hyperbolic conjugate points
# fig2d_sphere = conjugate.plot(['hyperbolic_with_loops_north', 'hyperbolic_with_loops_south',
# 'hyperbolic_without_loops'], linewidth=1.0, view=PLANE,
# color=color_conjugate, figure=fig2d_sphere, linestyle='solid')
# for fig3d_sphere in figs3d_sphere:
# conjugate.plot(['hyperbolic_with_loops_north', 'hyperbolic_with_loops_south',
# 'hyperbolic_without_loops'], linewidth=1.0, view=SPHERE,
# color=color_conjugate, figure=fig3d_sphere, linestyle='solid')
# add abnormal geodesics until the cusp point defined by the boundary of the domain of strong current
α1a, α2a = α_abnormals(r0, v)
#
print('α1a = ', α1a)
print('α2a = ', α2a)
# plot
fig2d_sphere = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=LW_CUT,
color=color_abnormal_2D, figure=fig2d_sphere,
zorder=geometry2d.plottings.z_order_splitting)
fig2d_sphere = geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=PLANE, linewidth=LW_CUT,
color=color_cut_locus, linestyle='dotted', figure=fig2d_sphere,
zorder=geometry2d.plottings.z_order_splitting)
for fig3d_sphere in figs3d_sphere:
geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=LW_CUT,
color=color_abnormal_3D, figure=fig3d_sphere,
zorder=geometry2d.plottings.z_order_splitting)
geodesic.plot(alphas=[α1a, α2a], tf=t_cusp_abnormal, view=SPHERE, linewidth=LW_CUT,
color=color_cut_locus, linestyle='dotted', figure=fig3d_sphere,
zorder=geometry2d.plottings.z_order_splitting)
α1a = 2.1639961029445254 α2a = 4.119189204235061
fig2d_sphere.savefig('kepler_spheres_2D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
fig2d_sphere
figs3d_sphere[0].savefig('kepler_spheres_3D.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
figs3d_sphere[0]
gap = 1e-3
delta = 0.025
rs = np.arange(gap, np.pi-gap, delta)
prs = np.arange(-1.0, 1.0, 1*delta)
RS, PRS = np.meshgrid(rs, prs)
PTS = np.zeros((len(prs), len(rs)))
def pthetafun(r, pr, pzero, i, j):
t = 0.0
q = np.array([r, 0.0]) # theta is a cyclic variable
def fun(ptheta):
p = np.array([pr, ptheta])
return prob.Hamiltonian(t, q, p) + pzero
# solve fun(ptheta) = 0
# for ptheta_guess, use a value close to PTS[i, j]: do not access PTS[i, j] directly or out of bounds
m, n = PTS.shape
if i == 0:
i_ = 1
elif i == m - 1:
i_ = m - 2
else:
i_ = i+1
if j == 0:
j_ = 1
elif j == n - 1:
j_ = n - 20
else:
j_ = j+1
pheta_guess = PTS[i_, j_]
#ptheta_sol = scipy.optimize.fsolve(fun, pheta_guess)
options = nt.nle.Options(Display='off')
sol = nt.nle.solve(fun, pheta_guess, options=options)
return sol.x
pzero = -1.0
for j, r in enumerate(rs):
for i, pr in enumerate(prs):
H = pthetafun(r, pr, pzero, i, j)
PTS[i,j] = H
fig, ax = plt.subplots()
CS = ax.contour(RS, PRS, PTS)
#ax.clabel(CS, inline=True, fontsize=10)
#ax.set_title('Simplest default with labels')
pass
delta = 0.01
prmax = 3.0
gap = 1e-3
rs = np.arange(0+gap, np.pi-gap, delta)
prs = np.arange(-prmax, prmax, delta/prmax)
RS, PRS = np.meshgrid(rs, prs)
def pthetafun(r, pr, pzero, sign):
m2_ = g2_Kepler(r)
v_ = v_Kepler
a = (1/m2_) - v_**2
b = -2*v*pzero
c = pr**2 - pzero**2
delta = b**2 - 4*a*c
if delta < 0.0:
return np.nan
else:
ptheta1 = (-b + np.sqrt(delta))/(2*a)
ptheta2 = (-b - np.sqrt(delta))/(2*a)
if -pzero - v_ * ptheta1 >= 0.0 and -pzero - v_ * ptheta2 >= 0.0:
if sign == 1:
m = max(ptheta1, ptheta2)
if m >= 0:
return m
else:
return np.nan
elif sign == -1:
m = min(ptheta1, ptheta2)
if m <= 0:
return m
else:
return np.nan
else:
raise geometry2d.errors.ArgumentValueError('sign must be equal to 1 or -1')
elif -pzero - v_ * ptheta1 >= 0.0:
m = ptheta1
if sign * m >= 0:
return m
else:
return np.nan
elif -pzero - v_ * ptheta2 >= 0.0:
m = ptheta2
if sign * m >= 0:
return m
else:
return np.nan
else:
return np.nan
# hyperbolic
pzero = -1.0
#
PTS_HYP_NEG = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
for i, pr in enumerate(prs):
pθ = pthetafun(r, pr, pzero, -1)
PTS_HYP_NEG[i,j] = pθ
PTS_HYP_POS = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
for i, pr in enumerate(prs):
pθ = pthetafun(r, pr, pzero, 1)
PTS_HYP_POS[i,j] = pθ
# elliptic
pzero = 1.0
PTS_ELL = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
for i, pr in enumerate(prs):
pθ = pthetafun(r, pr, pzero, -1)
PTS_ELL[i,j] = pθ
# abnormal
pzero = 0.0
PTS_ABN = np.zeros((len(prs), len(rs)))
for j, r in enumerate(rs):
for i, pr in enumerate(prs):
pθ = pthetafun(r, pr, pzero, -1)
PTS_ABN[i,j] = pθ
restart_calcul_level_lines = True
if prob.data.contains('level_lines') and not restart_calcul_level_lines:
level_lines = prob.data.get('level_lines')
PTS_HYP_NEG = np.array(level_lines['PTS_HYP_NEG'])
PTS_HYP_POS = np.array(level_lines['PTS_HYP_POS'])
PTS_ELL = np.array(level_lines['PTS_ELL'])
PTS_ABN = np.array(level_lines['PTS_ABN'])
print('level_lines loaded from data')
else:
level_lines = {'PTS_HYP_NEG': PTS_HYP_NEG.tolist(), 'PTS_HYP_POS': PTS_HYP_POS.tolist(),
'PTS_ELL': PTS_ELL.tolist(), 'PTS_ABN': PTS_ABN.tolist()}
prob.data.update({'level_lines': level_lines})
print('level_lines saved in data')
level_lines saved in data
# function to plot the domain of strong current between two angles r1 and r2
def geometry2d.utils.plot_2d_domain_level_lines(ax, r1, r2):
# plot the domain: a 2d surface
ylims = ax.get_ylim()
x = np.linspace(r1, r2, 100)
y1 = np.ones_like(x)*ylims[0]
y2 = np.ones_like(x)*ylims[1]
ax.fill_between(x, y1, y2, color=color_strong_domain_2d, alpha=0.5)
return ax
#
fig, ax = plt.subplots()
#
#CS = ax.contour(RS, PRS, PTS_HYP_NEG, [-4, -3, -2, -1, -0.25], colors=color_hyperbolic_2D, linewidths=1)
#ax.clabel(CS, inline=True, fontsize=10)
CS = ax.contour(RS, PRS, PTS_HYP_POS, [0.02, 0.25, 0.5, 0.75], colors=color_hyperbolic_2D, linewidths=1, linestyles='dashed')
#ax.clabel(CS, inline=True, fontsize=10)
# add geodesics
r0 = np.pi/2.0
θ0 = 0.0
tf = 6.12
pzero = -1.0
pr0s = np.linspace(1.25, 3.0, 4)
pθ0s = []
for pr0 in pr0s:
pθ0 = pthetafun(r0, pr0, pzero, -1)
pθ0s += [pθ0]
q0 = np.array([r0, θ0])
p0 = np.array([pr0, pθ0])
t0 = 0.0
N = 200
tspan = list(np.linspace(t0, tf, N+1))
qs, ps = prob.extremal(t0, q0, p0, tspan)
r = [q[0] for q in qs]
pr = [p[0] for p in ps]
ax.plot(r, pr, color=color_hyperbolic_2D, linewidth=1.0, linestyle='solid')
#
ax.set_xlim(0.0, np.pi)
ax.set_ylim(-3.0, 3.0)
geometry2d.utils.plot_2d_domain_level_lines(ax, r1_strong, r2_strong)
#
print('pθ0s = ', pθ0s)
fig.savefig('kepler_level_lines_hyperbolic.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
#
pass
pθ0s = [-0.3228913274337241, -1.1266351099358802, -1.9642005576625727, -2.8172904669025316]
#
fig, ax = plt.subplots()
#
CS = ax.contour(RS, PRS, PTS_ELL, [-20, -8, -6, -4, -3], colors=color_elliptic_2D, linewidths=1, linestyles='solid')
geometry2d.utils.plot_2d_domain_level_lines(ax, r1_strong, r2_strong)
if False:
# add geodesics
θ0 = 0.0
pr0 = 0.0
tf = 5
pzero = 1.0
r0s = np.linspace(1.1, 1.5, 4)
pθ0s = []
for r0 in r0s:
pθ0 = pthetafun(r0, pr0, pzero, -1)
pθ0s += [pθ0]
q0 = np.array([r0, θ0])
p0 = np.array([pr0, pθ0])
print('q0 = ', q0)
print('p0 = ', p0)
t0 = 0.0
N = 200
tspan = list(np.linspace(t0, tf, N+1))
qs, ps = prob.extremal(t0, q0, p0, tspan)
r = [q[0] for q in qs]
pr = [p[0] for p in ps]
ax.plot(r, pr, color=color_hyperbolic_2D, linewidth=1.0, linestyle='solid')
fig.savefig('kepler_level_lines_elliptic.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
#
pass
#
fig, ax = plt.subplots()
#
CS = ax.contour(RS, PRS, PTS_ABN, [-6, -4, -3, -2, -1, -0.25], colors='green', linewidths=1, linestyles='solid')
geometry2d.utils.plot_2d_domain_level_lines(ax, r1_strong, r2_strong)
if False:
# add geodesics
θ0 = 0.0
r0 = np.pi/2.0
tf = 5
pzero = 0.0
pr0s = np.linspace(0.1, 2.5, 4)
pθ0s = []
for pr0 in pr0s:
pθ0 = pthetafun(r0, pr0, pzero, -1)
pθ0s += [pθ0]
q0 = np.array([r0, θ0])
p0 = np.array([pr0, pθ0])
print('q0 = ', q0)
print('p0 = ', p0)
t0 = 0.0
N = 200
tspan = list(np.linspace(t0, tf, N+1))
qs, ps = prob.extremal(t0, q0, p0, tspan)
r = [q[0] for q in qs]
pr = [p[0] for p in ps]
ax.plot(r, pr, color=color_hyperbolic_2D, linewidth=1.0, linestyle='solid')
fig.savefig('kepler_level_lines_abnormal.pdf', dpi=300, transparent=True, format="pdf", bbox_inches='tight')
#
pass
# 3D plots
from matplotlib import cm
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
surf = ax.plot_surface(RS, PRS, PTS_HYP_POS, cmap=cm.coolwarm, linewidth=0, antialiased=False)
#ax.set_zlim(-20.0, 0.0)
min(PTS_HYP_NEG[1000,:])
-202.83241871800467